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We demonstrate a close analogy between a viscoelastic medium and an electrically con- 
ducting fluid containing a magnetic field. Specifically, the dynamics of the Oldroyd-B 
^ ' fluid in the limit of large Deborah number corresponds to that of a magnetohydrody- 

namic (MHD) fluid in the limit of large magnetic Reynolds number. As a definite example 
of this analogy, we compare the stability properties of differentially rotating viscoelas- 
. tic and MHD flows. We show that there is an instability of the Oldroyd-B fluid that is 

physically distinct from both the inertial and elastic instabilities described previously in 
i the literature, but is directly equivalent to the magnetorotational instability in MHD. It 

I occurs even when the specific angular momentum increases outwards, provided that the 

. angular velocity decreases outwards; it derives from the kinetic energy of the shear flow 

Oh' and does not depend on the curvature of the streamlines. However, we argue that the 

Q . elastic instability of viscoelastic Couette flow has no direct equivalent in MHD. 



1. Introduction 

1.1. Viscoelastic and magnetohydrodynamic fluids 
In his investigation of the viscosity of gases, Clerk Maxwell (1867) proposed that the 
stress in a fluid obeys an equation of the form 

, , d(stress) , . . . d(strain) _ , 

(stress) +T — (viscosity) x , (1.1) 

dt dt 

where t is the relaxation time. If the time scale of the straining motion is long compared 
to T, the second term on the left-hand side is negligible and the stress is proportional 
to the rate of strain. This Newtonian relation gives rise to viscous behaviour. However, 
if the time scale of the strain is short compared to r, the first term on the left-hand 
side is negligible and the stress is proportional to the strain itself. This Hookean relation 
gives rise to elastic behaviour. The reason for the elastic response is that the rapid 
strain prevents the configuration of the molecules from relaxing towards an equilibrium 
distribution, and the stress is therefore 'frozen in' to the fluid. 

Modern constitutive equations for viscoelastic fluids (Bird, Armstrong & Hassager 
1987a) are usually expressed in a covariant tensorial form based on the principles set out 
by Oldroyd (1950). His hquid B is one of the most widely used nonhnear models of a 
viscoelastic fluid, and provides a fair representation of a dilute solution of a polymer of 
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high molecular weight. It is based on Maxwell's equation (^]^), but cast in a form that 
satisfies the principle of material frame indifference. Moreover, it can also be derived from 
the kinetic theory of idealized extensible polymer molecules contained in a Newtonian 
solvent (Bird, Curtiss & Armstrong 1987&). The dimensionless number characterizing 
the ratio of the relaxation time to the time scale of the flow is the Deborah number (or 
Weissenberg number); when this is large, the polymeric stress is effectively 'frozen in' to 
the fluid. 

In an electrically conducting fluid, the magnetic field B affects the dynamics through 
the bulk Lorentz force (e.g. Roberts 1967). This can be represented in terms of the 
Maxwell electromagnetic stress tensor ,0 

M=^-^i, (1.2) 

the two parts of which correspond to a tension in the field lines and an isotropic magnetic 
pressure. It is well known that, in a perfectly conducting fluid, the magnetic field is 'frozen 
in' to the fluid, in the sense that magnetic field lines can be identified with material lines 
(Alfven 1950). Even in a fluid of finite conductivity, the magnetic field is effectively 
'frozen in' for motions of sufficiently short time scale, or sufficiently large length scale, 
corresponding to a large magnetic Reynolds number. It follows that the Maxwell stress 
is also 'frozen in' to the fluid in a certain sense. 

From a mathematical point of view, the 'freezing in' of a tensor field X(r,t) in a flow 
with velocity field M(r, t) can be expressed by the equation (e.g. Tur & Yanovsky 1993) 

r)X 

— + £^X = 0, (1.3) 
ot 

where £ is the Lie derivative. For a scalar field X, this gives the familiar expression 
dX/dt + u • 'VX = 0, meaning that the numerical value of X is conserved by every fluid 
element. For a (contravariant) vector field B it gives the induction equation of ideal, 
incompressible magnetohydrodynamics (MHD), which implies that the magnetic field is 
advected and stretched in the same way as infinitesimal line elements. For a second-rank 
tensor field it results in the 'upper-convected derivative' that appears in the governing 
equation of the Oldroyd-B fluid (Bird et al. 1987a). 

At a physical level, an analogy is to be seen between a polymer solution, containing 
extensible molecules that are advected and distorted by the flow and react on it through 
their tension, and an electrically conducting fluid, containing magnetic field lines that 
are also advected and distorted by the flow and react on it through their tension. 

It follows that there is a physical and mathematical similarity between the dynamics of 
viscoelastic and MHD fluids. We will show that a formal analogy can be drawn between 
the Oldroyd-B fluid in the limit of large Deborah number and an MHD fluid in the limit 
of large magnetic Reynolds number. In other words, in this limit, the Maxwell stress in 
MHD obeys the equation of a Maxwell fluid. 

1.2. Instabilities of differentially rotating fluids 
Differentially rotating flows are common in astrophysics and geophysics, and have been 
studied extensively in the laboratory. The simplest form of differential rotation occurs 
when the angular velocity depends only on the cylindrical radius, fl — 0(r), and Couette 
flow between differentially rotating cylinders provides an excellent model system for 
investigating the dynamics of such flows. According to Rayleigh (1916), instability occurs 

f Here is the permeability of free space. For non-relativistic flows, the electric fleld makes 
a negligible contribution to the stress. 
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in the absence of viscosity whenever the specific angular momentum |r^f2| decreases 
outwards. Most subsequent theoretical and experimental studies, starting with the classic 
work of Taylor (1923), have focused on the onset of Rayleigh's inertial instability in the 
presence of viscosity, and the interesting sequence of dynamical states that ensues. 

Numerous variants of Couette flow have also been considered, and among these the 
Couette flow of viscoelastic fluids has received much attention. Early work in the 1960s 
(e.g. Thomas & Walters 1964; Giesckus 1966) examined the effect of viscoelasticity on 
the onset of Rayleigh's inertial instability. However, one of the most important recent 
results is the theoretical and experimental demonstration by Larson, Shaqfeh & MuUer 

(1990) of a physically distinct instability in viscoelastic Couette flow. This is a purely 
elastic instability that occurs at sufHciently large Deborah number T|dn/dlnr|, even in 
the limit of negligible inertia, and irrespective of the sign of the angular momentum 
gradient. 

The influence of a magnetic field on the stability of the Couette flow of an electrically 
conducting fluid has also been investigated theoretically and (to a lesser extent) experi- 
mentally. Again, early work (described by Chandrasekhar 1961) focused on the effect of 
the magnetic field on the onset of inertial instability. More importantly, Velikhov (1959) 
and Chandrasekhar (1960) uncovered a physically distinct instability in magnetized Cou- 
ette flow. In the absence of viscosity and resistivity, and in the presence of a weak vertical 
magnetic fleld, this 'magnetorotational' instability occurs whenever the angular velocity 
decreases outwards, irrespective of the sign of the angular momentum gradient. 

The magnetorotational instability finds its most important applications in astrophysi- 
cal fluid dynamics, where magnetic fields are prevalent and the astronomical length scales 
allow for large magnetic Reynolds numbers. The stability of differentially rotating flows is 
of considerable interest in astrophysics, especially in connection with accretion discs (e.g. 
Pringle 1981). These are usually thin discs of gas in circular orbital motion around a star 
or black hole. The angular velocity decreases outwards according to Kepler's third law, 
fl cx r~^/^, and the Reynolds numbers are extremely high (e.g. lO'^'*). Observations indi- 
cate that angular momentum is transported outwards through accretion discs at a much 
greater rate than allowed by viscosity, and understanding the origin of this 'anomalous 
viscosity' has been a major goal of accretion disc research. 

The anomalous viscosity is usually attributed to turbulent transport. However, despite 
the very high Reynolds numbers, there is no convincing demonstration of any suitable 
hydrodynamic instability in circular Keplerian flow. Indeed, simple reasoning can be 
used to argue that hydrodynamic turbulence is unlikely to be self-sustaining in a flow 
that amply satisfles Rayleigh's stability criterion (Balbus & Hawley 1998). However, it 
has been demonstrated that MHD turbulence develops very readily in accretion discs, 
through the nonlinear development of the magnetorotational instability. Since the results 
of Velikhov (1959) and Chandrasekhar (1960) were rediscovered by Balbus & Hawley 

(1991) and their significance was appreciated, the magnetorotational instability has been 
analysed in considerable detail in the astrophysical literature. 



1.3. Properties of the magnetorotational instability 

The properties of the magnetorotational instability have been reviewed by Balbus & 
Hawley (1998), and we recall some of the important features here. Its simplest manifes- 
tation is in an incompressible, inviscid, perfectly conducting fluid having angular velocity 
Q{r) and containing a uniform magnetic fleld B parallel to the axis of rotation. An ap- 
proximate local dispersion relation can be obtained for normal modes having growth rate 
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s and wavevector k parallel to B, and has the form 

+ [AVLiVL - A) + 2cji] + uj\{uj\ - AQ.A) = 0, (1.4) 

where the quantity 

--if 

measures the differential rotation, and is known as Oort's first constant in the astrophys- 
ical literature. The Alfven frequency is wa = (a'op)""'^^^^ • -B, and the combination 

m{n- A) = ^^{r^n') (1.6) 

is the Rayleigh discriminant, or the square of the epicyclic oscillation frequency. When 
this is positive, unstable normal modes with s > can nevertheless be found provided 
that A^A > 0, as is the case in astrophysical discs. The maximal growth rate, s = \A\^ is 
achieved by a mode having lo\ = A{2Vl — A). 

More generally, and from a local perspective, rotation and shear in the correct rel- 
ative orientation are required, and a weak magnetic field of any geometry is sufficient 
to initiate the instability, provided the fiuid is sufficiently ionized. In the presence of 
significant dissipation, the ideal growth rate must compete with viscous and resistive 
damping, so that growth rates less than |^| are achieved, or the instability may be sup- 
pressed altogether. An unstable mode must always bend the field lines, having a non-zero 
Alfven frequency, and therefore the instability of a purely azimuthal (or toroidal) field is 
essentially non-axisymmetric. 

A simple explanation of the instability can be given in terms of two fluid elements, 
connected by magnetic field lines, that are initially in circular orbit at the same radius. 
The fiuid elements are then given angular momentum perturbations of opposite sign. The 
one receiving the positive perturbation moves to an orbit of larger radius and acquires 
a smaller angular velocity, lagging behind its partner. The tension of the magnetic field 
exerts a torque that pulls the lagging clement forwards, enhancing the initial perturbation 
and leading to instability. A mechanical analogue, consisting of two orbiting particles 
connected by a weak spring, also exhibits instability. 

1.4. Plan of the paper 

The main purpose of this paper is to draw attention to the physical and mathematical 
similarity between viscoelasticity and MHD. As an example of the application of this 
idea, we explore in some detail the relation between the instabilities of differentially 
rotating viscoelastic and MHD flows. As described in § |l.2| , previous investigations have 
uncovered instabilities of viscoelastic and MHD Couette flow that are physically distinct 
from Rayleigh's inertial instability. In the light of the analogy we describe, an obvious 
question is whether the elastic instability of Larson et al. (1990) is somehow related to 
the magnetorotational instability. We will argue that this is not the case, but will show 
that there is another instability of viscoelastic Couette flow that is the direct equivalent 
of the magnetorotational instability. 

The remainder of this paper is organized as follows. In § 2 we set out the basic equations 
governing incompressible viscoelastic and MHD flows and present the analogy between 
them. We then discuss, in § 3, the possible sources of instability on the basic of energy 
considerations. In § 4 we define a model system consisting of plane Couette flow in a 
rotating channel, equivalent to cylindrical Couette flow in the narrow-gap limit, and for- 
mulate eigenvalue problems for the normal modes of the Oldroyd-B and MHD fluids. 
Some numerical solutions are presented in § 5 to illustrate the expected similarity be- 
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tween the two systems. In § 6 the existence of locahzed growing solutions in the two 
systems, satisfying the same magnetorotational dispersion relation, is demonstrated by 
an asymptotic analysis. Finally, the results are summarized and discussed in § 7. 



2. Basic equations 

2.1. Viscoelastic fluid 

We first consider an incompressible viscoelastic fluid of uniform density p. The Oldroyd-B 
model is characterized by a solvent viscosity /i, a polymer viscosity /Xp and a relaxation 
time T. The velocity field U obeys the solenoidal condition, 

V-U = 0, (2.1) 

and the equation of motion 

— + U-Vu\=-V^ + V-T + iiV^U. (2.2) 

Here Vl/ = p + p$ is the modified pressure (<i> being the gravitational potential) and T 
is the Oldroyd-B stress, which is a sym metric tensor field of second rank satisfying the 
constitutive equation (cf. equation (1.1)) 



T - 



^ + U-VT- (VUf -T-T-VU 

ot 



- A^p [V(7 + (VJ7)T] , (2.3) 



where the superscript 'T' denotes the transpose of a second-rank tensor. 

2.2. MHD fluid 

We also consider an incompressible, electrically conducting fluid of uniform density p, 
viscosity p, and electrical conductivity a. The velocity field U obeys the solenoidal con- 
dition, 

V • J7 0, (2.4) 

and the equation of motion 

— + U-Vu\=-V^ + V-M + p\I^U. (2.5) 

Here 5* = p + /?$ is again the modified pressure and M is the Maxwell stress given in 
equation ( |l.2| ). The magnetic field obeys the solenoidal condition, 

V-J5 = 0, (2.6) 

and the induction equation, 

— + U -VB = B -VU + j]V^B, (2.7) 

which is derived from Maxwell's equations for the electromagnetic field and Ohm's law. 
Here t] = l/{poa) is the magnetic diffusivity. 

2.3. The formal analogy 

Instead of comparing the Oldroyd-B stress T and the Maxwell stress M directly, we take 
the polymeric and magnetic stress tensors to be 

Tp=T+^J, T,„ = — , (2.8) 

T Mo 
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which differ from T and M only by the addition of an isotropic part in each case. 
Specifically, Tp does not include the equilibrium isotropic pressure nkT of the poly- 
mer molecules, and Tm does not include the isotropic pressure /2^q of the magnetic 
field. In an incompressible fluid, such terms can be taken care of by writing the two 
equations of motion in the form 

or r \ 

U-VU]= - V*p,,„ + V • Tp,,„ + ^iV^U, (2.9) 



dt 



where 



^'p = «' + ^, 4,^^^, + ^ (2.10) 

are redefined modified pressures. 

According to the Oldroyd-B constitutive equation (2.3) and the induction equation 
( ^.7[ ), Tp and satisfy the equations 



dt 



^^ + i7.VTp-(Vt/)T. Tp - Tp.VC/ = -^(rp-^i), (2.11) 



O -7- 

— f + c/ • vr„, - {vuf • r„, - t„, • vt/ = ^ [bs/^b + {v^b)b] . (2.12) 

dt Mo 
In the limits r — > oo and 77 — > 0, corresponding to large Deborah number and large mag- 
netic Reynolds number respectively, the right-hand sides of these equations are negligible. 
The polymeric and magnetic stresses then satisfy identical equations, involving the same 
upper-convected derivative, and they appear identically in the equation of motion of the 
fluid. Therefore the formal analogy can be expressed symbolically as 

lim (Oldroyd-B fluid) = lim(MHD fluid). (2.13) 

r — ^00 rj — >0 

We note that is a positive semi-definite tensor having one non-negative eigenvalue 
and two zero eigenvalues. Joseph (1990) has shown t hat T p also retains a positive definite 



character when it evolves according to equation ( 2.11 ). This is required on physical 
grounds, because in the derivation of the Oldroyd-B constitutive equation from kinetic 
theory, Tp cx (dd), where d is the separation of the ends of a pol ymer molecule, and 



the angle brackets denote an average (Bird et al. 1987&). Equation (|2.11 ) shows that Tp 



attempts to return to isotropy on the relaxation time, but in a shear flow at large Deborah 
number, this tendency is overcome and one eigenvalue of Tp does indeed dominate, as 
required by the MHD analogy. 

As we have shown, the induction equation of ideal MHD provides an equation for the 
magnetic stress tensor, which is comparable to the constitutive equation of the Oldroyd- 
B fluid. One might ask whether the constitutive equation can be reduced to something 
resembling an induction equation. This is indeed so: at any instant we may express the 
positive definite tensor Tp in terms of three vector fields Bi, 

1 3 

Tp^—Y^B^B,. (2.14) 



Equation ( 2.11 ) is recovered if the fields evolve according to the induction- like equations 
dB, 



dt 



U-VB, = B,.VU-^[B,-i^Q,>j, (2.15) 
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1 ^ 

-Y,{B,Q, + Q,B,)^1. (2.16) 

i=l 

To find the fields Qj, let B be the ma trix whose columns are (Bi, B2, -B3), and similarly 
for Q. In matrix notation, equation ( 2.16| ) reads 

\{BQ^ + QST) ^ J (2.17) 

and is satisfied when 

Q^{1+A)C, (2.18) 
where A is an arbitrary antisymmetric matrix and C is the inverse of . This means 
that 

= C, + OxC„ (2.19) 
where Jl is an arbitrary vector field. The vector fields Ci are just the reciprocal vectors 
to {Bi], e.g. 

The non-uniqueness of the fields reflects the fact that the representation ( ^.14[ ) is 
partially redundant: we are expressing a tensor field with six independent components 
in terms of three vector fields each having three independent components. It is therefore 
permissible to impose three constraints on the fields B^, and then will no longer be 
arbitrary. For example, it may be convenient to require that the fields Bi be solenoidal. 
Now equation (2.15) implies that 

I + J7 . v) (V . B.) = (V • B, - ^ V . q) , (2.21) 

and so the solenoidal property is preserved if is chosen such that 

= V • g, = V • (Cj + f2 X C), i = l,2,3. (2.22) 
Choosing the fields Bi to be solenoidal also ensures that 



1 ^ 

V- Tp = — VS, -VB,, (2.23) 
Mo — 



for direct comparability with the Lorentz force. 



3. Energetics and instability 

3.1. MHD fluid 

Some insight into the possible instabilities of viscoelastic and MHD flows can be obtained 
on the basis of energy considerations. Instabilities typically release energy stored in the 
basic state and use this to allow a perturbation to grow in time. This restricts the class 
of flows that can exhibit instability, and limits the growth rates that can be achieved. 

We start by considering the case of the MHD fluid, which is more straightforward. 
Starting from the equations of § 2.2 it is possible to derive an energy equation of the 
form 

- + V.F^-A (3.1) 
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where 



is the energy density, 



E^\PU^ + ^ (3.2) 



F = (E + -—{B- U)B -fiUx{V xU)- x {V x B) (3.3) 
Mo Mo 



is the energy flux, and 

D^/ilVxJJp + ^IVxBp (3.4) 
Mo 

is the dissipation rate. 

Consider a perturbative solution of the equations of § 2.2 in which upper-case symbols 
([/, B, 4') denote the basic state (not necessarily steady) and lower-case symbols (it, b, ^fj) 
denote the Eulerian perturbations. Using the linearized equations, it is then possible to 
derive the energy-like equation 

I- (hu^ + ^) + V • F' = (-puu + —) ■.VU+{bxu)-J 

at V 2^0/ V Mo/ 



~^l\V xu\^ - xb\^ (3.5) 

Mo 

governing the perturbations, where F' is a certain flux and J = Mo x B is the current 
density in the basic state. The quantity differentiated with respect to time is the part 
of the energy density at second order in the perturbation amplitude that must grow in 
any instability. Provided that the instability is local, so that it does not depend on a 
particular choice of boundary conditions, the term V • F' cannot play an essential role 
in this equation, because it will vanish on integration over the volume of the fluid in the 
case of periodic boundary conditions or, in many cases, physical boundary conditions. 
Therefore any local instability must derive its energy either from the kinetic energy of 
the basic flow, through the term involving VfJ, or from the magnetic energy, through the 
term involving J . For kinetic energy to be released, there must be a velocity gradient, 
because a uniform flow can be eliminated by a Galilean transformation and therefore 
cannot be a source of instability. A potential magnetic field (J — 0) also cannot be 
a source of instability, as it minimizes the magnetic energy in a region subject to the 
magnetic flux through its boundary being prescribed (e.g. Priest 1982). 

In the case of a potential magnetic field it is possible to place an upper bound on the 
growth rate of any local instability. Let 

S = i [VJ7 + ( VC/)T] (3.6) 

be the rate-of-strain tensor of the basic fiow, and (Ai, A2, A3) its eigenvalues. Its quadratic 
form satisfies the inequalities 

S \ XX 

min(Ai, A2,A3) ^ — '— — s$ max(Ai, A2, A3) (3.7) 



and therefore 



^max(|Ai|,|A2|,|A3|). (3.8) 



It follows from equation (3.5) that the largest possible growth rate of any local instability 
is the largest eigenvalue, in absolute value, of the rate-of-strain tensor of the basic flow. 
When the magnetic held is not potential the growth rate can be increased by at most 

(/.o/m)^/^|J|/2. 
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3.2. Viscoelastic fluid 

By working with the representation ( ^.14 ) of the polymeric stress in terms of three vector 
fields Bi it is possible to derive a similar energy-like equation 



,2 




- V • F" = i-puu + y ^1 : Vt/ + y (6, xu)-J, 

+ ±J2i^.b,)u.B,-^,\Vxuf--J2^ + -^i2^^-Q^ (3.9) 

governing linear perturbations from any basic state, where F" is a certain flux and 
Ji = /ip X Bi by analogy with MHD. If we constrain the representation such that 
the fields are solenoidal, then the term involving V • bi vanishes. The argument proceeds 
almost as before, with gradients in the basic fiow or the basic stress providing potential 
sources of energy for the disturbance. The final term, involving 6j • q^, is a third possible 
source of energy, but the dependence suggests that the effect of this term may be 
expected to be small in the limit of large Deborah number. 



4. Plane Couette flow in a rotating channel 

As a minimal model of a differentially rotating fiow, we consider a linear shear flow 
(plane Couette fiow) in a rotating channel. This is equivalent to cylindrical Couette flow 
in the limit of a narrow gap, if the angular velocities of the two cylinders are not widely 
disparate. All effects of curvature are then neglected. 

4.1. Basic state and boundary conditions 

We adopt Cartesian coordinates (x, y, z) in a frame of reference rotating with uniform 
angular velocity . The only change required to the equations in the rotating frame 
is the inclusion of the Coriolis force. The centrifugal force, which is derivable from a 
potential, can be absorbed into the modified pressure. The rotation of the frame does 
not affect any of the other equations. The equation of motion therefore becomes 

p(^^ + U -VU + 2ne^ xU^ ^ ■■■ . (4.1) 

We consider fiow in the channel < x < d between a stationary plane boundary 
a: = and a moving plane boundary x = d with velocity —2 Ad By. The non-slip and 
impermeable boundary conditions 

U = Q at a; = 0, U = -2Adey at x^d (4.2) 

apply. The basic flow is the plane Couette fiow, 

U = -2Axey. (4.3) 

A modified pressure quadratic in x is required to balance the Coriolis force. 

When this model is taken as a local representation of a differentially rotating flow 
with angular velocity ri(r), the shear parameter A is to be interpreted as Oort's first 
constant. When the Rayleigh discriminant 4ri(ri — A) is positive, the flow of an inviscid, 
unmagnetized flow is linearly stable to axisymmetric perturbations. For a Keplerian flow 
in which Q. oc r~'^/^, we have A/^ = 3/4. (In the rheological literature, the shear rate 
\2A\ would usually be called 7.) 
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In the case of the viscoelastic fluid, the non-zero stress components associated with 
the basic flow are 



P' 



(4.4) 



which provide the steady solution of equation (2.3). The polymeric stress defined in 
equation (|2.8|) is 



1 



-De 







-De 2De^ + l 











1 



(4.5) 



where De = 2At is the Deborah number, and this can be represented in the form (2.14) 
using the three solenoidal fields 



Bi 



1/2 



-1 

De ± {De^ + 1)^/2 




(4.6) 



Note that, for large De, the field -Bi is much greater than the other two and corresponds 
to a uniform magnetic field almost exactly in the y-direction. 

In the case of the MHD fiuid, we suppose that a uniform magnetic field B = By By 
is imposed. We also suppose the boundaries to be perfectly conducting, so that the 
additional boundary conditions 



dBy ^ dB, 

dx dx 







(4.7) 



apply at a; = and x = d. 

We note that the magnetic stress tensor in the MHD fluid resembles the polymeric 
stress tensor in the viscoelastic fluid if De is large and we identify 



Mo 



8AV^p. 



(4.8) 



The energy considerations of § ^ are not affected by the rotation of the frame of 
reference, because the Coriolis force does no work on the fluid. The eigenvalues of the 
rate-of-strain tensor are {A, — A, 0). As the magnetic field is uniform, the maximal growth 
rate of any local instability, at least in the MHD case, is \A\. Incidentally, this proves 
the conjecture of Balbus & Hawley (1992) that the magnetorotational instability, with a 
suitably chosen wavevector and in the absence of dissipation, achieves the largest possible 
growth rate of any local shear instability. (In an inviscid, unmagnetized fiuid, the largest 
possible growth rate of Rayleigh's inertial instability is AQ,{A — il). This is always less 
than or equal to \A\, with equality in the case A — 2f2.) 



4.2. Dimensionless groups 

We introduce the kinematic viscosities v — ii/p and ^'p — fJ-p/p. The four dimensionless 
parameters of the viscoelastic system are the Rossby number, the Reynolds number, the 
Deborah number and the viscosity ratio, defined by 

Ro = ^, Re^——, De = 2AT, S = (4.9) 

respectively. For reference, the Taylor number is Ta = Re^Ro^^{l — Ro^^). 

The four dimensionless parameters of the MHD system are the Rossby number, the 
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Reynolds number, the magnetic Reynolds number and the Chandrasekhar number, de- 
fined by 

A 2Ad^ 2Ad'^ BU^ , ^ 

Ro^-, Re = , Rm^ , Q=— — , 4.10 

respectively. The identification (^|^) corresponds to 

g • (4.11) 

We can therefore quote the magnetic field strength in terms of an effective Deborah 
number for the MHD system, to make a direct comparison easier. 

We are interested in comparing the behaviour of the viscoelastic and MHD systems in 
the limit of large De and large Rm. This limit could be approached in many different 
ways, but we choose to do this as one might in an ideal experiment, by keeping p, u, Up, 
T, 7] and d fixed while increasing f2, A and By together. This means that Ro and S are 
fixed while Re cx Rm oc De and Q oc De^. The elasticity, De/ Re, is fixed in this process. 

4.3. Linear perturbations 

We now consider small deviations from the above state, such that the Eulerian pertur- 
bation of velocity, say, is 

Re [u{x) exp(si + \kyy + ifczz)] , (4-12) 

where s is the growth rate (in general complex) and ky and kz are real wavenumber 
components. Differentiation of the perturbations with respect to x will be denoted by a 
prime, and we define k"^ = ky + k'^. 

4.4. Viscoelastic fluid 
The perturbations of the viscoelastic fluid satisfy the equations 



u'^ + ikyUy + \kzUz — 0, 
p [{s - 2\Axky)ux - 2Q.Uy\ = + t'^^ + ikyt^y + ifcz^xz + /^(w" - k'^u^), 

P [(S - 2\Axky)Uy + 2{n - A)U,JC\ = -ikyip + t'^y + \kytyy + ik^tyz + p{Uy - k^ Uy) , 

p{s - 2iAxky)uz = -ikz^J + t'^^_ + ikytyz + ik^tzz + ^^{Uz ~ ^^"z), 
txx + [(s ~ 2iAxky)txx ~ 2iTxykyUx] = 2ppU^, 
txy + [{s 2iAxky)txy -\- 2Atxx ^xyij^x ^ ^^y^y) ~ ^TyykyUx] — PpiUy + ikyU 



txz + [(s ~ 2\Axky)txz ~ ^TxykyUz] — Ppiuz -\- ikzUx), 
iyy ^ [("5 2iAxky^tyy -\- AAtxy 2Txy^y ^^-^yy^y^y] — 2ippkyUy , 

tyz + [(s 2iAxky)tyz + 2Atxz ^ TxyUz '^T^yy^y'^z\ — ^fJ-pi^yUz -|- kzUy), 

tzz +t[{s ~ 2iAxky)tzz] = 2ippkzUz- 
These constitute a sixth-order system of linear ODEs to be solved for the eigenvalue s. 
The dependent variables may be taken as {ip, Ux,Uy,Uy,Uz,u'z). The boundary conditions 
Ux — Uy — Uz — apply at a; = 0, d. 

These equations must be solved numerically in general. However, it is instructive to 
analyse further the case of unsheared or 'axisymmetric' modes (ky = 0), which correspond 



4.13) 

4.14) 
4.15) 
4.16) 
4.17) 

4.18) 
4.19) 

4.20) 

4.21) 

4.22) 
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to axisymmetric modes in cylindrical geometry. In this case, the equations have constant 
coefficients and can be combined into a single equation for Ux^ 

[qs + [qv + Dux + 4f^(f^ - A)q^klux - AnATi^pk^Dux = 0, (4.23) 

where q = I + ts and D is the operator 

D.-i^+e (4.24) 

This equation may be investigated analytically in an approximate way by considering 
solutions of a simple trigonometric form Ux oc shi{kxx), although these cannot satisfy all 
six physical boundary conditions. (In § ^ below we compute the global solutions of this 
equation numerically.) The local dispersion relation corresponding to these solutions, 

[qs + {qv + Vp){kl + kl)f{kl + kl)+^n{^-A)q^kl-A^ATVpkl{kl + kl) = 0, (4.25) 

is a quartic equation for s with real coefficients. It can be shown that the principle of the 
exchange of stabilities holds: instability first sets in at a stationary bifurcation (s = 0), 
which occurs when the constant term passes through zero, i.e. when 

[v + Vpf{kl + klf + AVL{VL - A)kl - AVlATVpkl{kl + kl) = 0. (4.26) 

Suppose that Rayleigh's criterion for stability, 451 (fJ — A) > 0, is satisfied, and that 
AVlA > 0. When t is increased from zero to a sufficiently large value, a bifurcation occurs 
and axisymmetric instability ensues. To understand this we note that, when v = 0, and 
in the limit r ^ with fc^ ^ fc^, the dispersion relation ( 4.25| ) becomes identical 

to the ideal magnetorotational dispersion relation ( |l.4| ) for a vertical magnetic field and 
vertical wavevector, provided that we identify ^ ^Q^-p/r. This is precisely what 
is suggested by the field of equation ( |4.6| ). Although the principal analogy is with 
a uniform magnetic field in the y-direction, such a field provides no restoring force to 
axisymmetric perturbations and we see instead the effect of the much weaker field B^. 
Therefore the axisymmetric viscoelastic instability, which we verify numerically in § ^ 
below, can be understood as being analogous to a magnetorotational instability deriving 
from the weak vertical field. 

4.5. MHD fluid 
The perturbations of the MHD fluid satisfy the equations 

u'^ + ikyUy + ikzUz = 0, (4.27) 

p \{s - 2\Axky)ux - 2nuy] = -V'm + ^^J-o^kyBybx + n{Ux - k^Ux)., (4.28) 

p [[s — 2\Axky)uy + 2(il — A)ux] = — ifcyV'm + '^t^o^^y^v^y + l^i^y ~ k'^Uy), (4.29) 

p{s — 2\Axky)uz = — ifcz?Am + ifJ-o^kyBybz + m("z ^ k'^u^), (4.30) 

b'x + ikyby + ikzbz ^ 0, (4.31) 

(s - 2iAxky)bx = ikyByUx + r]{b'^ - k'^bx), (4.32) 

(s - 2iAxky)by + 2Abx = ikyByUy + r]{b'^ - k'^by), (4.33) 

(s - 2iAxky)bz = ikyByUz + viK - k'^^z)- (4.34) 
These constitute a tenth-order system of linear ODEs to be solved for the eigenvalue 
s. The dependent variables may be taken as {ipai,Ux,Uy,u'y,UzTu'z,by^ b' bz,b'z). The 
boundary conditions Ux = Uy = Uz = by = b'z = apply at x = 0, d. To eliminate bx 
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from the problem, differentiate equation (4.31) to find then substitute into equation 
(PH) to find 



(s - 2\Axky + rik^)bx — ikyByUx - T]{ikyb' + ikzb'J. 



(4.35) 



Therefore b^ is determined algebraically in terms of the dependent variables, and can be 
substituted where needed. It automatically satisfies the boundary condition b^ = at 
X = 0, d. A difficulty would arise if the quantity s — 2iAxky + rjk"^ were to vanish at any 
point. As this is a complex function, it is 'unlikely' that both real and imaginary parts 
would vanish simultaneously. In any case, it could vanish only for a decaying mode, and 
such modes are of no interest here. 



Although the linearized equations (4.27)-(4.34) for the MHD fluid appear quite dif- 
ferent from those of the viscoelastic fluid, equations (4.13)-( 4.22| ), they can be seen to 



correspond in the limits r — > oo, ^ if we identify Tyy ^ By/ fig, t^y ^ Byb^/na 



Bybz/fj-o, 



while Txyj txx^ ixz^ tzz 



0. 



tyy -(-^ 2Byby/ Ho aud tyz 

In the special case of axisymmetric modes {ky = 0) the velocity and magnetic pertur- 
bations are decoupled. The magnetic perturbation always decays if ?7 > 0. The remaining 
equations can be combined into a single equation for Ux-, 



(s + vDfT)Ux + 4fi(fi - A)k^Ux = 0. 



(4.36) 



Stability is assured if Rayleigh's criterion, — A) > 0, is satisfied. 



5. Numerical investigation 

We solve the eigenvalue problems defined in §§ 4^ and 4^ for non- axisymmetric modes 
numerically by the shooting method. The arbitrary normalization "0(0) = 1 is adopted, 
and the equations are integrated from a; = to x = 1. For the viscoelastic system, the 
boundary conditions at a; = 1 impose three conditions on the three unknown quantities 
s, My(0) and u^(0). Newton-Raphson iteration is applied to converge on a solution. For 
the MHD system, shooting in is required. 

In the absence of viscosity and resistivity, the MHD problem becomes identical to 
the 'Cartesian model' studied by Ogilvie & Pringle (1996) in their investigation of the 
magnetorotational instability in the presence of an azimuthal (or toroidal) magnetic field. 
We recall some results of that analysis: (i) the instability requires a non-zero azimuthal 
wavenumber ky so that the magnetic field lines are bent by the perturbation; (ii) as kz is 
increased, unstable modes emerge from the continuous spectrum of Alfven waves and the 
eigenvalues approach limit points; (iii) the largest growth rates are attained in the limit 
kz — > oo, when the normal modes are localized near a boundary (although, as we show 
below, solutions also exist that grow rapidly but transiently in the interior of the fiuid); 
(iv) the maximal growth rate. A, is attained for an Alfven frequency lj^ = (15/16) '^/^Jl 
in the Keplerian case A/Vl. = 3/4. 

For numerical purposes it is convenient to adopt d and il"^ as units of length and time. 
We adopt Ro = 3/4, which is stable according to Rayleigh's criterion and is suggested 
by astrophysical applications, and take v — v-p = rj for simplicity. 

In the presence of dissipation, all modes decay in the limit kz ^ oo. Therefore we 
restrict attention to a moderate value, kzd = tt, at which the growth rates are appreciable 
(but not optimal, and always less than A). We select the optimal Alfven frequency, as 



described in § 1.3, by choosing kyd = 1 and De/Re — 5/24. The most unstable mode is 
one with the fewest nodes in its eigenfunction. The variation of its growth rate with De 
is shown in figure 1. As De increases, the eigenvalues of the mode in the viscoelastic and 
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Figure 1. Variation of the growth rate of the first unstable non-axisymmetric mode (at ky = 1) 
with the Deborah number, for the viscoelastic fluid (solid line) and the MHD fluid (dashed line). 
Note that the Reynolds number is Re = (24/5)-De. 



MHD systems converge. The eigenfunctions are also in close agreement at De = 150, as 
shown in figure 2. 

We have also solved numerically for axisymmetric unstable modes in the viscoelastic 
problem, as anticipated in § 4.4. The growth rates of the most unstable modes are shown 
in figure 3. They are smaller than for the non-axisymmetric modes, consistent with the 
idea that the axisymmetric instability is analogous to a magnetorotational instability 
deriving from the weak vertical field Bj,. 



6. Asymptotic analysis 

We now present an asymptotic analysis demonstrating the existence of localized, non- 
axisymmetric growing solutions of the perturbation equations for both systems, consistent 
at leading order with the magnetorotational dispersion relation ( |l.4| ). We are interested 
again in the limit of large Deborah number and large magnetic Reynolds number. 

6.1. MHD fluid 

We consider a solution of the perturbation equations for the MHD fluid, localized in a 
layer near an arbitrary point x — xq. Let 



X = xq + eX, 



(6.1) 
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Figure 2. Eigonfunctions of a non-axisymmctric unstable mode at De = 150, for the viscoelastic 
fluid (left) and the MHD fluid (right). Real and imaginary parts are shown with solid and dashed 
lines, respectively. 



where e <C 1 is an ordering parameter and X = 0(1) within the layer of interest. We 
introduce the scalings 

k, = e-^/%, y = e^i), V = e\ (6.2) 

implying that the vertical wavelength is even shorter than the width of the layer, and 
that dissipation has only a weak effect on the solution. The solution will have an expo- 
nential time-dependence at leading order, but we relax the assumption of a normal mode 
and allow the solution to evolve freely on a long time scale captured by the slow time 
coordinate T = et. This is achieved through the replacement 



(6.3) 
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0.12 




Figure 3. Variation of the growth rates of axisymmetric unstable modes {ky = 0) with the 
Deborah number, for the viscoelastic fluid. The eigenfunctions of modes appearing successively 
as De is increased have increasing numbers of nodes. 



in the perturbation equations, and we also replace 



_d_ ^1 d 

within the layer. A consistent expansion scheme for the perturbations is of the form 



(6.4) 





= UM{X,T) + eu,i{X,T)- 




Uy 


= UyQ{X, T) + eUyi{X, T) ^ 






= ei/2 [u,,(X,T) + 0{e)], 






= e^ [MX,T) + 0{e)], 




bx 


= bM{X,T) + ebxi{X, T)4 


'0{e% 


by 


^byo{X,T) + ebyi{X,T)^ 




bz 


= £1/2 [bzo{X,T) + 0{e)]. 





(6.5) 



From equations (i.2S), (4.29), (4.32) and (4.33) at leading order we obtain the algebraic 
system 

p{suxo - 2^Uyo) = ifig^kyBybxo, 

p [sUyQ + 2{Q - A)uxo] = ifJ-Q^kyBybyO, 

sbxo = ikyByUxo, 
sbyo + 2AbM = ikyByUyo, (6.6) 
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where s — sq — 2iAkyXo. These may be combined into the single equation 

{s"^ + J2 -A) + 2ujI] + luKloI - mA)} = O, (6.7) 

where = kyBy/ ^op. This has a non-trivial solution, 

UM = F{X,T), (6.8) 

if and only if s satisfies the magnetorotational dispersion relation ( |l.4| ). Wc then d educe 
Uyo, &xO an d b„o in terms of F, and also u^o, V'o a-nd ^^o from equations (4.27), (4.30) 

and either ( 4.31 ) or ( 4.34 ) at le ad ing o rder. 

From equations (|4.2S| ), (|4.29|) , ( [4.32| ) and ( [4.33|) at order e we similarly obtain 



{s^ + [m{n -A) + 2u;l] + ujUljI - iflA)} u^i = R, 



(6.9) 



where the right-hand side R depends on F and its derivatives. Given that s has been 
chosen to satisfy the dispersion relation, the solvability condition for this equation is 
i? = 0, which results in an evolutionary equation for F, 



dF d'^F 

_.„_ + ,i«-c)K 

This is a modified diffusion equation containing constant coefficients 

{s+uir 



2/c? 



s|s^ + 2c^2 s2+^2(^2 ^4^2)]- 



2 j4/co/ 



(6.10) 



(6.11) 



(6.12) 



(6.13) 



s* + 2ijis^+u;iiiji + 4n^) 

When the conditions for instability are met, s is real and positive and therefore a, b and c 
are real and a and c are positive. A particular solution of equation ( 6.10| ), corresponding 
to an initial condition F{X, 0) = and valid for T > in the absence of boundaries, 

is the Green function 



F = (47rar)"^/^exp 



12acT2 - GiabXT'^ + 3X^ 
12ar 



(6.14) 



as can be obtained by Fourier-transform methods. The Green function decays as T cx) 
for any fixed X , or as \X\ —^ oo for any fixed T. 

It follows that localized solutions exist that grow exponentially at leading order, follow- 
ing the magnetorotational dispersion relation. The envelope of the solution evolves more 
slowly in time but ultimately decays superexponentially, so that the instability grows for 
many e-folding times before the development of very short length scales leads to decay. 
If we insisted on having a normal- mode solution, equation (6.10) would become an Airy 
equation in X. It can be shown (Ogilvie 1997) that localized solutions of this type do 
exist, but only near the boundaries of the fiuid. 

The reason for considering disturbances that are localized in x is that it provides a 
convenient method of demonstrating the existence of growing solutions without resorting 
to numerical analysis. Provided that the localization scale 6x is long compared to the 
vertical wavelength and to the characteristic dissipative scales, the growth rate is insen- 
sitive to Sx. Terquem & Papaloizou (1996) have shown that, in the limit of ideal MHD, 
a continuous spectrum of infinitely localized growing disturbances exists. 
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6.2. Viscoelastic fluid 

A very similar analysis can be carried out for the viscoelastic fluid. The additional re- 
quirements are the scaling t = e^^f and the expansions 

^vv ~ iyyoiX, T) + etyyi{X, T) + O(e^), 
ty.=e'^^ [ty,o{X,T) + 0{e)], 

t,.=0{e'). (6.15) 
Otherwise the analysis is so similar that we do not repeat it in detail. Equations 



and (6.9) are obtained exactly as before, provided that we identify = kyTyy/p in the 



dispersion relation. Exactly the same evolutionary equation (|6.10| ) is also obtained with 
the sole exception that the term involving fj does not appear in the coefficient c. 

Therefore this method also establishes the correspondence between the viscoelastic and 
MHD fluids in the limit of large De and large Rm, and demonstrates the existence of the 
magnetorotational instability in the Oldroyd-B fluid. 



7. Discussion 

We have demonstrated a close analogy between a viscoelastic medium and an electri- 
cally conducting fluid containing a magnetic field. Both an Oldroyd-B fluid, in the limit 
of large Deborah number, and a magnetohydrodynamic fluid, in the limit of large mag- 
netic Reynolds number, feature a stress tensor that is nearly 'frozen in' to the fluid in a 
precise mathematical sense. As a definite example of this analogy, we have examined a 
local model of a differentially rotating fluid, consisting of plane Couette flow in a rotating 
channel. The stress tensor in the case of a viscoelastic fluid resembles the Maxwell stress 
corresponding to a magnetic field aligned with the flow. 

Our analysis demonstrates that there is a detailed correspondence between instabili- 
ties in the two systems. We have identified a direct equivalent of the magnetorotational 
instability in the viscoelastic fluid. It exists when the angular velocity and relative vortic- 
ity are antiparallel (or when the angular velocity decreases outwards) and the maximal 
growth rate is equal to the shear parameter, or Oort constant, A. It is distinguished most 
clearly from Rayleigh's inertial instability by the fact that it occurs even when the spe- 
cific angular momentum increases outwards. It is also distinct from the elastic instability 
described by Larson et al. (1990), which depends on the curved geometry of Couette flow 
and exists in the elastic limit, De/ Re oo. 

We have also found an axisymmetric viscoelastic instability that can be understood 
as being analogous to the magnetorotational instability of a vertical magnetic field. This 
reflects the fact that the polymeric stress tensor Tp can be decomposed into three effective 
magnetic fields, one of which is a uniform field almost aligned with the flow and another 
of which is a uniform vertical field. Although the vertical field is much weaker in the limit 
of large -De, it provides the dominant restoring force for axisymmetric disturbances. 

The instability discussed by Larson et al. (1990) does not appear in our analysis be- 
cause we have neglected the curvature of the streamlines. Although Larson et al. (1990) 
considered the limit of a narrow gap, the characteristic growth rates they obtained are 
smaller than the growth rates we have discussed, by a factor of order e^/^, where e is the 
ratio of the gap width to the radius. Being purely elastic in nature, the instability of Lar- 
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son et al. (1990) must derive its energy from the elastic energy stored in the flow, rather 
than the shear energy which is the source for inertial and magnetorotational instabilities. 

It is natural to enquire whether there is an MHD equivalent of the instability of Larson 
et al. (1990), in which energy is derived from the magnetic field. In cylindrical Couette 
flow at large De the dominant stress component has the form T^^ oc r~^, which can be 
identified with an azimuthal (or toroidal) magnetic field oc r~^. Typically, toroidal 
pinch configurations arc unstable to modes that rely on the curvature of the magnetic 
field lines and derive energy from the magnetic configuration. The most dangerous are 
the m = 1 'kink' modes and m = 'sausage' modes. However, in the absence of fluid 
motion the profile oc is sufficiently steep to be stable to all perturbations (Tayler 
1973). If the profile of T^^ happened to be less steep, for example T^^ oc r~^, it is likely 
that there would be a viscoelastic equivalent of the kink instability. We therefore conclude 
that the instability of Larson et al. (1990) is not directly related to an MHD instability, 
but relies on inherently viscoelastic effects not captured by our analogy. This conclusion 
is supported by an examination of the physical explanation that Larson et al. give for 
their instability. 

Since the work of Larson et al. (1990) there have been a number of related theoretical 
and experimental studies of viscoelastic Couette flow, some of which have examined the 
effects of inertia and non-axisymmetry (e.g. Avgousti & Boris 1993; Steinberg & Grois- 
man 1998; Baumert & MuUer 1999). Some of these might have been expected to reveal 
the analogue of the magnetorotational instability. However, it appears that the cases 
usually investigated are those in which cither the outer or inner cylinder is stationary, or 
the inner cylinder rotates at twice the angular velocity of the outer cylinder with only 
a narrow gap between the two. This means that, whenever the analogue of the magne- 
torotational instability might have occurred, the system is unstable to Raylcigh's inertial 
instability. In order to separate the two effects, it would be valuable to examine cases 
in which the angular velocity decreases outwards but the speciflc angular momentum 
increases. Interestingly, the magnetorotational instability has never been demonstrated 
in laboratory experiments. Although there is currently much effort towards this goal 
(e.g. Goodman & Ji 2002), the technical requirements are considerable and the system is 
constrained by the very small magnetic Prandtl numbers of liquid metals. A viscoelastic 
magnetorotational experiment might prove to be less demanding and easier to visualize. 

We conclude with some further perspectives on the analogy between viscoelastic and 
MHD flows. 

The analogy is asymptotic in nature and therefore not perfect. Viscoelastic and MHD 
flows deviate from simple stress freezing in different ways: the viscoelastic stress relaxes, 
while the magnetic fleld diffuses. The classes of exactly steady solutions of the two systems 
are therefore different, because after a sufficiently long time cither relaxation or diffusion 
will have its effect. For example, while there is only one solution for viscoelastic Couette 
flow, magnetized Couette flow can be set up with either vertical or azimuthal current- 
free magnetic fields. In this sense, the analogy is more applicable to dynamical, time- 
dependent situations than to steady fiows. 

Renardy (1997) has analysed the large-De limit of steady, two-dimensional flows of 
the upper-convcctcd Maxwell fluid (obtained by setting the solvent viscosity fi of the 
Oldroyd-B fluid to zero). Noting that the stress T has one dominant eigenvalue in this 
limit, he writes T = puu, where p and u are a flctitious density and a flctitious velocity 
field. When inertial forces arc negligible, V • T must balance the pressure gradient and 
p and u are then found to satisfy the steady Euler equations. Through our analogy, a 
connection can be seen here with the work of Moffatt (1985), which makes use of the 
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analogy between steady Euler flows and magnetostatic equilibria in which the Lorentz 
force balances the pressure gradient. 

There is a close connection between the polymeric stress Tp at moderate De and the 
mean stress tensor {BB)/^o of a disordered magnetic field, as occurs in MHD turbu- 
lence. Under these conditions, Tp and (BB) /jio may have three positive eigenvalues of 
comparable magnitude. Indeed, Ogilvie (2001) has suggested the use of Maxwellian vis- 
coelastic models, with Deborah numbers of order unity, for MHD turbulence in accretion 
discs. 

Besides Couette flow, another problem that has received much attention is the stabil- 
ity of a planar jet or shear layer with respect to two-dimensional disturbances. Azaiez 
& Homsy (1994) and Rallison & Hinch (1995) examined the equivalent problem for a 
viscoelastic fluid, noting the potentially stabilizing influence of a polymer additive. Tak- 
ing a limit in which Re and De tend to infinity while maintaining a finite ratio, they 
derived an elastic equivalent of Raylcigh's stability equation. An analogy exists between 
this problem and that of the stability of a similar flow of an ideal, electrically conducting 
fluid with a magnetic field parallel to the flow, a problem studied since the 1950s. In a 
recent study, Hughes & Tobias (2001) derived a magnetic Rayleigh equation (their equa- 
tion (3.5)) that is exactly equivalent to the elastic Rayleigh equation given by Rallison 
& Hinch (1995, p. 314). 

In this paper we have drawn attention to a useful analogy between viscoelastic and 
MHD flows, and have discussed the relation bc;twecn instabilities of differential rotation 
in the two systems. We anticipate that much further use can be made of this analogy. 

We are grateful to John Hinch, John Rallison and the referees of this paper for their 
constructive comments. G.I. O. acknowledges the support of the Royal Society through 
a University Research Fellowship. 
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